% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for simulations w/ external predictors
% Called by: simsExtPred.m
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function simOut = simsExtPredFcn(s,data2,data,p,sp)

    rng(12345+s);

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Simulate data
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [simfwd,simPredX] = simulateData(data2,data,p,sp,2);

    simyld = cumsum(simfwd,2)./cumsum(ones(size(simfwd)),2);
    simdret = zeros(size(simfwd));
    simdret(1:p.H,:) = NaN;
    simdret(p.H+1:end,2:end) = simfwd(1:end-p.H,2:end) - simfwd(p.H+1:end,1:end-1);
    simxret = cumsum(simdret,2);

    simX = simfwd(:,p.factorFwdMats);
    simX1 = [ones(sp.T,1),simX];
    
    betaHatAll = NaN(numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    betaSEHatAll = NaN(numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    betaSEHatAll_LLSW_EWC = NaN(numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    betaSEHatAll_LLSW_NW = NaN(numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));

    for i = 1:numel(p.yldMats)
        tmp = nwest(simxret(p.H+1:end,i),[simX1(1:end-p.H,:) simPredX(1:end-p.H,:)],p.H);
        betaHatAll(i,:) = tmp.beta;
        betaSEHatAll(i,:) = tmp.beta./tmp.tstat;
        tmp = nwest(simxret(p.H+1:end,i),[simX1(1:end-p.H,:) simPredX(1:end-p.H,:)],sp.S_LLSW_NW);
        betaSEHatAll_LLSW_NW(i,:) = tmp.beta./tmp.tstat;
        tmp = ewc(simxret(p.H+1:end,i),[simX1(1:end-p.H,:) simPredX(1:end-p.H,:)],sp.nu_LLSW_EWC);
        betaSEHatAll_LLSW_EWC(i,:) = tmp.betahat./tmp.that;
    end    

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BH Bootstrap
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    [simfacload,simevals] = OrderedPCA(simyld, 'cov');
    Y1 = simyld*simfacload(:,1:3);
    simehat = simyld - Y1*simfacload(:,1:3)';
    sims2hat = trace(simehat'*simehat)/numel(simehat);

    bsOutBH = bootstrapVAR(Y1,simPredX,sp.B,p.BH_B_BC,p.BH_B_BC,1,1);
    
    %!% Initialize Matrices
    bsBetaHatAll_CG = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    bsBetaSEHatAll_CG = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    bsBetaHatAll_BH = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    bsBetaSEHatAll_BH = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));   
    bsBetaHatAll_FullOracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));
    bsBetaSEHatAll_FullOracle = NaN(sp.B,numel(p.yldMats),size(data2.X1,2)+size(data.pred,2));

    for b = 1:sp.B
    
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % BH Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        bsYldBH = bsOutBH.bsY1(:,:,b)*simfacload(:,1:3)' + sqrt(sims2hat)*randn(sp.T,numel(p.yldMats));
        bsFwdBH = [bsYldBH(:,1) bsYldBH(:,2:end).*(ones(sp.T,1)*p.yldMats(2:end))/data.period-bsYldBH(:,1:end-1).*(ones(sp.T,1)*p.yldMats(1:end-1))/data.period];
    
        bsDretBH = zeros(size(bsYldBH));
        bsDretBH(1:p.H,:) = NaN;
        bsDretBH(p.H+1:end,2:end) = bsFwdBH(1:end-p.H,2:end) - bsFwdBH(p.H+1:end,1:end-1);
        bsXretBH = cumsum(bsDretBH,2);

        bsX_BH = bsFwdBH(:,p.factorFwdMats);
        bsX1_BH = [ones(sp.T,1),bsX_BH];
        
        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretBH(p.H+1:end,i),[bsX1_BH(1:end-p.H,:), bsOutBH.bsY2(1:end-p.H,:,b)],p.H);
            bsBetaHatAll_BH(b,i,:) = tmp.beta;
            bsBetaSEHatAll_BH(b,i,:) = tmp.beta./tmp.tstat;
        end

        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Oracle Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        [bsFwdFullOracle,bsPredFullOracle] = simulateData(data2,data,p,sp,2);

        bsDretFullOracle = zeros(size(bsFwdFullOracle));
        bsDretFullOracle(1:p.H,:) = NaN;
        bsDretFullOracle(p.H+1:end,2:end) = bsFwdFullOracle(1:end-p.H,2:end) - bsFwdFullOracle(p.H+1:end,1:end-1);
        bsXretFullOracle = cumsum(bsDretFullOracle,2);

        bsX_FullOracle = bsFwdFullOracle(:,p.factorFwdMats);
        bsX1_FullOracle = [ones(sp.T,1),bsX_FullOracle];
        
        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretFullOracle(p.H+1:end,i),[bsX1_FullOracle(1:end-p.H,:), bsPredFullOracle(1:end-p.H,:)],p.H);
            bsBetaHatAll_FullOracle(b,i,:) = tmp.beta;
            bsBetaSEHatAll_FullOracle(b,i,:) = tmp.beta./tmp.tstat;
        end       
        
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % CG Bootstrap
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        simdretapprox = simfwd(1:end-1,2:end) - simfwd(2:end,1:end-1);
    
        varJoint = var1BC([simfwd(:,end) simPredX],p.CG_B_BC);
        Zfwd = varJoint.VBC(2:end,1);
        Zpred = varJoint.VBC(2:end,2:end);
        Zinit = [Zfwd, simdretapprox, Zpred];
        Z = [Zinit; Zinit(1:sp.M-1,:)];
        idxRand = randi(sp.T-sp.M,[ceil(sp.T/sp.M),1]);
        idxBootstrap = kron(idxRand,ones(sp.M,1)) + kron(ones(ceil(sp.T/sp.M),1),(0:sp.M-1)');
        idxBootstrap = idxBootstrap(1:sp.T-1);
        bsZ = Z(idxBootstrap,:);
        bsDret = [NaN(1,size(simfwd,2)); NaN(size(simfwd,1)-1,1) bsZ(:,2:numel(p.yldMats))];        
        %
        bsV = bsZ(:,[1,end-p.L+1:end]);    
        bsV = bsV - ones(sp.T-1,1)*mean(bsV);        
        tmpJoint = NaN(sp.T,p.L+1);
        tmpJoint(1,:) = [simfwd(1,end) simPredX(1,:)];
        for t = 2:sp.T
            tmpJoint(t,:) = [1 tmpJoint(t-1,:)]*varJoint.PsiBC + bsV(t-1,:);
        end
    
        bsFwdN = tmpJoint(:,1);
        bsPred_CG = tmpJoint(:,end-p.L+1:end);
        bsFwdCG = buildForwards(bsDret, simfwd(1,:), bsFwdN);        
        %{
        % Check buildForwards
        bsFwd2 = NaN(sp.T,numel(p.yldMats));
        bsFwd2(:,end) = filter(1 ,[1 -varFwd.PhiBC],...
               [0; varFwd.muBC+(bsZ(:, 1)-mean(bsZ(:, 1)))],simfwd(1,end));
        bsDret2 = bsZ(:,2:numel(p.yldMats));
        bsFwd2(1,:) = simfwd(1,:);
        for n = numel(p.yldMats)-1:-1:1
            for t = 2:sp.T
                bsFwd2(t,n) = bsFwd2(t-1,n+1) - bsDret2(t-1,n);        
            end
        end
        %}
        
        bsYldCG = cumsum(bsFwdCG,2)./cumsum(ones(size(bsFwdCG)),2);
        bsDretCG = zeros(size(bsFwdCG));
        bsDretCG(1:p.H,:) = NaN;
        bsDretCG(p.H+1:end,2:end) = bsFwdCG(1:end-p.H,2:end) - bsFwdCG(p.H+1:end,1:end-1);
        bsXretCG = cumsum(bsDretCG,2);

        bsX_CG = bsFwdCG(:,p.factorFwdMats);
        bsX1_CG = [ones(sp.T,1),bsX_CG];                

        for i = 1:numel(p.yldMats)
            tmp = nwest(bsXretCG(p.H+1:end,i),[bsX1_CG(1:end-p.H,:) bsPred_CG(1:end-p.H,:)],p.H);
            bsBetaHatAll_CG(b,i,:) = tmp.beta;
            bsBetaSEHatAll_CG(b,i,:) = tmp.beta./tmp.tstat;
        end
    end

simOut.betaHatAll = betaHatAll;
simOut.betaSEHatAll = betaSEHatAll;
simOut.bsBetaHatAll_CG = bsBetaHatAll_CG;
simOut.bsBetaSEHatAll_CG = bsBetaSEHatAll_CG;
simOut.bsBetaHatAll_BH = bsBetaHatAll_BH;
simOut.bsBetaSEHatAll_BH = bsBetaSEHatAll_BH;
simOut.bsBetaHatAll_FullOracle = bsBetaHatAll_FullOracle;
simOut.bsBetaSEHatAll_FullOracle = bsBetaSEHatAll_FullOracle;
simOut.betaSEHatAll_LLSW_EWC = betaSEHatAll_LLSW_EWC;
simOut.betaSEHatAll_LLSW_NW = betaSEHatAll_LLSW_NW;
simOut.simEigvals = simevals;
